########################### PERENNIAL AUSTERITY: DO LPS MAKE A DIFFERENCE? ###################

# Fiscal consolidations are implemented in the form of budgetary plans. These plans usually encompass not only the present budgetary
# period, but instead span over multiple years. For this reason, politicians have an incentive to delay the retrenchment of popular
# welfare state expenditures into the future, so as to shift the impact of austerity to future period and, perhaps, future incumbents.
# If we estimate the budgetary impact of an aggregated spending-side consolidation as defined in the paper, the initial impact
# of the spending adjustment will disproportionately concern other expenditure items and not social spending. For this reason, 
# dynamic long-run multipliers that are derived from the initial shocks as estimated in the first two periods are likely biased when 
# it comes to the true, nonmonotonic impact of austerity over the long run, here multiple (up to 8) years. 
#
# Leveraging the unique dataset provided by Alesina et al. (2019), I demonstrate three key findings: First, I show descriptively that 
# the retrenchment of social spending, on average, occurs at a later stage of the multiyear adjustment plan compared with public
# investment, which is almost exclusively cut in the initial period of the plan. Second, I demonstrate for Sweden the United Kingdom, 
# countries where the median social spending cut in the dataset occurs in the third year of the adjustment plan, that the canonical ARDL 
# with lag order (1,1) shows no significant short- or long-run impact of fiscal adjustment on social spending. In contrast, LPs evidence
# a persistent contraction in insurance-related social spending over time - exactly what the strategic timing of cutbacks argument
# predicts. Third, I show that autoregressive models with an extended lag order perform much better and provide very similar results
# as standard local projections. 

library(readxl)
library(tidyverse)

##### A) TOTAL ADJUSTMENT SHARES ######

alesinaplan <- read_xlsx("./Raw Files/Alesina Plan Incidence.xlsx")

# select all welfare state items, including eduaction spending: 
#  health, pensions, unemployment, other social spending, family and children

alesinaplan$wstate <- ifelse(alesinaplan$Category=="EDU"|alesinaplan$Category=="HLT"|alesinaplan$Category=="PENS"
                             |alesinaplan$Category=="UNEM"|alesinaplan$Category=="OTHSS"|alesinaplan$Category=="FCPO", 1, 0)

alesinaplan$invest <- ifelse(alesinaplan$Category=="INV"|alesinaplan$Category=="RD", 1, 0)

alesinaplan <- alesinaplan %>% dplyr::filter(wstate==1|invest==1)

alesinaplan$planid <- paste0(alesinaplan$Country, sep="_", alesinaplan$Year, sep="_", alesinaplan$measure)

alesinaplan <- alesinaplan %>% dplyr::select(planid,wstate,impact_t,impact_t1,impact_t2,impact_t3,impact_t4,impact_t5)

alesinaplan <- reshape2::melt(alesinaplan, id=c("planid","wstate"))

alesinaplan <- alesinaplan %>% dplyr::filter(!is.na(wstate) & value>=0)

alesinaplan$time <- NA

alesinaplan$time[alesinaplan$variable=="impact_t"] <- "0"
alesinaplan$time[alesinaplan$variable=="impact_t1"] <- "1"
alesinaplan$time[alesinaplan$variable=="impact_t2"] <- "2"
alesinaplan$time[alesinaplan$variable=="impact_t3"] <- "3+"
alesinaplan$time[alesinaplan$variable=="impact_t4"] <- "3+"
alesinaplan$time[alesinaplan$variable=="impact_t5"] <- "3+"

alesinaplan <- alesinaplan %>% group_by(wstate,time) %>%
  dplyr::mutate(value = sum(value, na.rm=T)) %>% ungroup() %>%
  distinct(wstate,time, .keep_all = T) %>%
  group_by(wstate) %>% dplyr::mutate(impact_sum = sum(value, na.rm=T)) %>%
  ungroup() %>% dplyr::mutate(impactsh = value/impact_sum)

p1 <- ggplot(data=alesinaplan, aes(x=time, y=impactsh, fill=factor(wstate))) +
  stat_summary(fun="sum", geom="bar", position = position_dodge(0.9)) +
  scale_fill_manual(values = c("0"="steelblue3","1"="red3"),
                    labels = c("0"="Public Investment", "1"="Welfare State")) +
  scale_y_continuous(expand = c(0,0)) +
  labs(fill="") + xlab("Year of Adjustment Plan") +
  ylab("Share of Total Adjustment Impact") +
  theme_pubclean() +
  theme(text = element_text(family="Times New Roman"),
        axis.line = element_line(color="grey", linewidth = 0.6),
        axis.ticks = element_line(color="grey"), legend.position = "right") 

ggsave(p1, filename="./Main Figures/adjustment_timing_sh.png", dpi=800, width=6, height=3.5)


#### B) BINNED YEARS TOTAL ADJUSTMENT #####

alesinaplan <- read_xlsx("./Raw Files/Alesina Plan Incidence.xlsx")

alesinaplan$wstate <- ifelse(alesinaplan$Category=="EDU"|alesinaplan$Category=="HLT"|alesinaplan$Category=="PENS"
                             |alesinaplan$Category=="UNEM"|alesinaplan$Category=="OTHSS"|alesinaplan$Category=="FCPO", 1, 0)

alesinaplan$invest <- ifelse(alesinaplan$Category=="INV"|alesinaplan$Category=="RD", 1, 0)

alesinaplan <- alesinaplan %>% dplyr::filter(wstate==1|invest==1)

alesinaplan$planid <- paste0(alesinaplan$Country, sep="_", alesinaplan$Year, sep="_", alesinaplan$measure)

alesinaplan <- alesinaplan %>% dplyr::select(planid,wstate,impact_t,impact_t1,impact_t2,impact_t3,impact_t4,impact_t5)

alesinaplan <- reshape2::melt(alesinaplan, id=c("planid","wstate"))

alesinaplan <- alesinaplan %>% dplyr::filter(!is.na(wstate) & value>=0) %>%
  group_by(planid) %>%
  dplyr::mutate(impact_sum = sum(value, na.rm=T)) %>%
  ungroup() %>% dplyr::mutate(impactsh = value/impact_sum)

alesinaplan$time <- NA

alesinaplan$time[alesinaplan$variable=="impact_t"] <- "0"
alesinaplan$time[alesinaplan$variable=="impact_t1"] <- "1"
alesinaplan$time[alesinaplan$variable=="impact_t2"] <- "2"
alesinaplan$time[alesinaplan$variable=="impact_t3"] <- "3+"
alesinaplan$time[alesinaplan$variable=="impact_t4"] <- "3+"
alesinaplan$time[alesinaplan$variable=="impact_t5"] <- "3+"

p1 <- ggplot(data=alesinaplan, aes(x=time, y=value, fill=factor(wstate))) +
  stat_summary(fun="sum", geom="bar", position = position_dodge(0.9)) +
  scale_fill_manual(values = c("0"="steelblue3","1"="red3"),
                    labels = c("0"="Public Investment", "1"="Welfare State")) +
  labs(fill="") + xlab("Year of Adjustment Plan") + ylab("Total Adjustment Impact") +
  theme_pubclean() + theme(text = element_text(family="Times New Roman"),
                           axis.line = element_line(color="grey", linewidth = 0.6),
                           axis.ticks = element_line(color="grey"), legend.position = "top") 

ggsave(p1, filename="adjustment_timing_sum.JPEG", dpi=500, width=5.5, height=3.5)


######### C) TOTAL ADJUSTMENT IMPACT RAW YEARS #########

alesinaplan <- read_xlsx("./Raw Files/Alesina Plan Incidence.xlsx")

alesinaplan$wstate <- ifelse(alesinaplan$Category=="EDU"|alesinaplan$Category=="HLT"|alesinaplan$Category=="PENS"
                             |alesinaplan$Category=="UNEM"|alesinaplan$Category=="OTHSS"|alesinaplan$Category=="FCPO", 1, 0)

alesinaplan$invest <- ifelse(alesinaplan$Category=="INV"|alesinaplan$Category=="RD", 1, 0)

alesinaplan <- alesinaplan %>% dplyr::filter(wstate==1|invest==1)

alesinaplan$planid <- paste0(alesinaplan$Country, sep="_", alesinaplan$Year, sep="_", alesinaplan$measure)

alesinaplan <- alesinaplan %>% dplyr::select(planid,wstate,impact_t,impact_t1,impact_t2,impact_t3,impact_t4,impact_t5)

alesinaplan <- reshape2::melt(alesinaplan, id=c("planid","wstate"))

alesinaplan <- alesinaplan %>% dplyr::filter(!is.na(wstate) & value>=0) %>%
  group_by(planid) %>%
  dplyr::mutate(impact_sum = sum(value, na.rm=T)) %>%
  ungroup() %>% dplyr::mutate(impactsh = value/impact_sum)

alesinaplan$time <- NA

alesinaplan$time[alesinaplan$variable=="impact_t"] <- 0
alesinaplan$time[alesinaplan$variable=="impact_t1"] <- 1
alesinaplan$time[alesinaplan$variable=="impact_t2"] <- 2
alesinaplan$time[alesinaplan$variable=="impact_t3"] <- 3
alesinaplan$time[alesinaplan$variable=="impact_t4"] <- 4
alesinaplan$time[alesinaplan$variable=="impact_t5"] <- 5

p1 <- ggplot(data=alesinaplan, aes(x=time, y=value, fill=factor(wstate))) +
  stat_summary(fun="sum", geom="bar", position = position_dodge(0.9)) +
  scale_fill_manual(values = c("0"="steelblue3","1"="red3"),
                    labels = c("0"="Public Investment", "1"="Welfare State\nSpending")) +
  labs(fill="") +
  xlab("Year of Adjustment Plan") +
  ylab("Total Adjustment Impact; % of GDP") +
  scale_x_continuous(breaks = c(0,1,2,3,4,5)) +
  theme_pubclean() +
  theme(text = element_text(family="Times New Roman"),
        axis.line = element_line(color="grey",
                                 linewidth = 0.6),
        axis.ticks = element_line(color="grey"),
        legend.position = "right") 

ggsave(p1, filename="adjustment_timing_sum.JPEG", dpi=500, width=6, height=3.5)



#### 2) THE CASES OF SWEDEN AND THE UNITED KINGDOM ######

countryplanner <- function(cntry){
  df <- read_xlsx("./Raw Files/Alesina Plan Incidence.xlsx")
  
  df$wstate <- ifelse(df$Category=="EDU"|df$Category=="HLT"|df$Category=="PENS"
                      |df$Category=="UNEM"|df$Category=="OTHSS"|df$Category=="FCPO", 1, 0)
  
  df$invest <- ifelse(df$Category=="INV"|df$Category=="RD", 1, 0)
  
  df <- df %>% dplyr::filter(Country==cntry) %>% dplyr::filter(wstate==1|invest==1)
  
  df$planid <- paste0(df$Country, sep="_", df$Year, sep="_", df$measure)
  
  df <- df %>% dplyr::select(planid,wstate,impact_t,impact_t1,impact_t2,impact_t3,impact_t4,impact_t5)
  
  df <- reshape2::melt(df, id=c("planid","wstate"))
  
  df <- df %>% dplyr::filter(!is.na(wstate) & value>=0) %>% group_by(planid) %>%
    dplyr::mutate(impact_sum = sum(value, na.rm=T)) %>%
    ungroup() %>% dplyr::mutate(impactsh = value/impact_sum)
  
  df$time <- NA
  
  df$time[df$variable=="impact_t"] <- 0
  df$time[df$variable=="impact_t1"] <- 1
  df$time[df$variable=="impact_t2"] <- 2
  df$time[df$variable=="impact_t3"] <- 3
  df$time[df$variable=="impact_t4"] <- 4
  df$time[df$variable=="impact_t5"] <- 5
  
  return(df)
}


### united kingdom ####

uk <- countryplanner("GBP")

p1 <- ggplot(data=uk, aes(x=time, y=value, fill=factor(wstate))) +
  stat_summary(fun="sum", geom="bar", position = position_dodge(0.9)) +
  scale_fill_manual(values = c("0"="steelblue3","1"="red3"),
                    labels = c("0"="Public Investment", "1"="Welfare State")) +
  labs(fill="") + xlab("Year of Adjustment Plan") + ylab("Total Adjustment Impact; % of GDP") +
  scale_x_continuous(breaks = c(0,1,2,3,4,5)) + theme_pubclean() +
  theme(text = element_text(family="Times New Roman"),
        axis.line = element_line(color="grey", linewidth = 0.6),
        axis.ticks = element_line(color="grey"), legend.position = "bottom") 

ggsave(p1, filename="uk_adjustment_timing_sum.JPEG", dpi=500, width=5, height=3.5)


### sweden ####

se <- countryplanner("SWE")

p2 <- ggplot(data=se, aes(x=time, y=value, fill=factor(wstate))) +
  stat_summary(fun="sum", geom="bar", position = position_dodge(0.9)) +
  scale_fill_manual(values = c("0"="steelblue3","1"="red3"),
                    labels = c("0"="Public Investment", "1"="Welfare State")) +
  labs(fill="") + xlab("Year of Adjustment Plan") + ylab("Total Adjustment Impact; % of GDP") +
  scale_x_continuous(breaks = c(0,1,2,3,4,5)) + theme_pubclean() +
  theme(text = element_text(family="Times New Roman"),
        axis.line = element_line(color="grey", linewidth = 0.6),
        axis.ticks = element_line(color="grey"), legend.position = "bottom") 

ggsave(p2, filename="se_adjustment_timing_sum.JPEG", dpi=500, width=5, height=3.5)



####################### REGRESSION ESTIMATES: ARDL & LOCAL PROJECTIONS #######################################

######### United Kingdom ##########

#### ARDL #####

macro_uk <- macro_lp %>% dplyr::filter(country=="United Kingdom")

macro_uk <- macro_uk %>%
  dplyr::mutate(ins_d0 = insurance_spending - dplyr::lag(insurance_spending, 1),
                ins_d1 = dplyr::lead(insurance_spending, 1) - dplyr::lag(insurance_spending, 1),
                ins_d2 = dplyr::lead(insurance_spending, 2) - dplyr::lag(insurance_spending, 1),
                ins_d3 = dplyr::lead(insurance_spending, 3) - dplyr::lag(insurance_spending, 1),
                ins_d4 = dplyr::lead(insurance_spending, 4) - dplyr::lag(insurance_spending, 1),
                ins_d5 = dplyr::lead(insurance_spending, 5) - dplyr::lag(insurance_spending, 1),
                inv_d0 = invest_spending - dplyr::lag(invest_spending, 1),
                ins_d0_l1 = dplyr::lag(ins_d0, 1),
                realgdpgr_l1 = dplyr::lag(realgdpgr, 1),
                realgdpgr_l2 = dplyr::lag(realgdpgr, 2),
                SPENDING_l1 = dplyr::lag(SPENDING, 1),
                unemp_l1 = dplyr::lag(unemp, 1),
                unemp_l2 = dplyr::lag(unemp, 2),
                insurance_spending_l1 = dplyr::lag(insurance_spending, 1))

m_ardl <- lm(insurance_spending ~ SPENDING + insurance_spending_l1 + SPENDING_l1 + realgdpgr + unemp + realgdpgr_l1 + unemp_l1 + realgdpgr_l2 + unemp_l2, data = macro_uk)

summary(m_ardl)

coef(m_ardl)

theta1 <- m_ardl$coefficients[3]
phi0 <- m_ardl$coefficients[2]
phi1 <- m_ardl$coefficients[4]

est1 <- phi0
est2 <- phi0 + theta1*phi0 + phi1
est3 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1)
est4 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2)
est5 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2) + (theta1^4*phi0 + phi1*theta1^3)
est6 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2) + (theta1^4*phi0 + phi1*theta1^3) + (theta1^5*phi0 + phi1*theta1^4)
estlr <- (phi0+phi1)/(1-m_ardl$coefficients[3])

se1 <- msm::deltamethod(~x2, mean = coef(m_ardl), cov = vcovHC(m_ardl))
se2 <- msm::deltamethod(~x2 + (x3*x2+x4), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se3 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se4 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3) + (x3^3*x2+x4*x3^2), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se5 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3) + (x3^3*x2+x4*x3^2) + (x3^4*x2+x4*x3^3), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se6 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3) + (x3^3*x2+x4*x3^2) + (x3^4*x2+x4*x3^3) + (x3^5*x2+x4*x3^4), mean = coef(m_ardl), cov = vcovHC(m_ardl))

Horizon <- c(1,2,3,4,5,6)
Estimate <- c(est1,est2,est3,est4,est5,est6)
se <- c(se1,se2,se3,se4,se5,se6)

impdat1 <- cbind(Horizon,Estimate,se)
impdat1 <- as.data.frame(impdat1)

impdat1$Estimate1 <- impdat1$Estimate
impdat1$Estimate2 <- impdat1$Estimate

impdat1$Estimate1[impdat1$Horizon>2] <- NA 
impdat1$Estimate2[impdat1$Horizon<=2] <- NA 

p1 <- ggplot(data=impdat1, aes(x=Horizon, y=Estimate1)) + geom_point(color="steelblue3") + geom_line(aes(y=Estimate), color="steelblue3") + geom_ribbon(aes(ymin=Estimate-se*1.645, ymax=Estimate+se*1.645), fill="steelblue3", alpha=0.2) + geom_hline(yintercept=0, linetype=2, linewidth=0.8, alpha=0.4) + ylab("Estimate") + scale_x_continuous(breaks = c(1,2,3,4,5,6), limits = c(1,6)) + scale_y_continuous(limits = c(-9,7)) + theme_pubclean() + theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "bottom")


#### LP ####

m0 <- lm(ins_d0 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_uk)
m1 <- lm(ins_d1 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_uk)
m2 <- lm(ins_d2 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_uk)
m3 <- lm(ins_d3 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_uk)
m4 <- lm(ins_d4 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_uk)
m5 <- lm(ins_d5 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_uk)

stargazer(m0,m1,m2,m3,m4,m5,type="text")

vcov0 <- coeftest(m0, vcov. = vcovHAC(m0, type="HC1"))
vcov1 <- coeftest(m1, vcov. = vcovHAC(m1, type="HC1"))
vcov2 <- coeftest(m2, vcov. = vcovHAC(m2, type="HC1"))
vcov3 <- coeftest(m3, vcov. = vcovHAC(m3, type="HC1"))
vcov4 <- coeftest(m4, vcov. = vcovHAC(m4, type="HC1"))
vcov5 <- coeftest(m5, vcov. = vcovHAC(m5, type="HC1"))

vcov0 <- vcov0[2,1:2]
vcov1 <- vcov1[2,1:2]
vcov2 <- vcov2[2,1:2]
vcov3 <- vcov3[2,1:2]
vcov4 <- vcov4[2,1:2]
vcov5 <- vcov5[2,1:2]

uklpdf <- as.data.frame(rbind(vcov0,vcov1,vcov2,vcov3,vcov4,vcov5))

uklpdf <- uklpdf %>% dplyr::rename(SE = `Std. Error`)

uklpdf$Horizon <- c(1:6)

p2 <- ggplot(data=uklpdf, aes(x=Horizon, y=Estimate)) + geom_point(color="red3") + geom_line(aes(y=Estimate), color="red3") + geom_ribbon(aes(ymin=Estimate-SE*1.645, ymax=Estimate+SE*1.645), fill="red3", alpha=0.2) + geom_hline(yintercept=0, linetype=2, linewidth=0.8, alpha=0.4) + scale_y_continuous(limits = c(-5,3)) + scale_x_continuous(breaks = c(1,2,3,4,5,6)) + ylab("LP Cumulative Multiplier") + theme_pubclean() + theme(text = element_text(family="Times New Roman"), axis.line = element_line(color="grey", linewidth = 0.6), axis.ticks = element_line(color="grey"), legend.position = "bottom")



###### SWEDEN ######


macro_se <- macro_lp %>% dplyr::filter(country=="Sweden")

macro_se <- macro_se %>%
  dplyr::mutate(ins_d0 = insurance_spending - dplyr::lag(insurance_spending, 1),
                ins_d1 = dplyr::lead(insurance_spending, 1) - dplyr::lag(insurance_spending, 1),
                ins_d2 = dplyr::lead(insurance_spending, 2) - dplyr::lag(insurance_spending, 1),
                ins_d3 = dplyr::lead(insurance_spending, 3) - dplyr::lag(insurance_spending, 1),
                ins_d4 = dplyr::lead(insurance_spending, 4) - dplyr::lag(insurance_spending, 1),
                ins_d5 = dplyr::lead(insurance_spending, 5) - dplyr::lag(insurance_spending, 1),
                ins_d6 = dplyr::lead(insurance_spending, 6) - dplyr::lag(insurance_spending, 1),
                ins_d7 = dplyr::lead(insurance_spending, 7) - dplyr::lag(insurance_spending, 1),
                inv_d0 = invest_spending - dplyr::lag(invest_spending, 1),
                ins_d0_l1 = dplyr::lag(ins_d0, 1),
                realgdpgr_l1 = dplyr::lag(realgdpgr, 1),
                realgdpgr_l2 = dplyr::lag(realgdpgr, 2),
                SPENDING_l1 = dplyr::lag(SPENDING, 1),
                unemp_l1 = dplyr::lag(unemp, 1),
                unemp_l2 = dplyr::lag(unemp, 2),
                insurance_spending_l1 = dplyr::lag(insurance_spending, 1))

m_ardl <- lm(insurance_spending ~ SPENDING + insurance_spending_l1 + SPENDING_l1 + realgdpgr + unemp + realgdpgr_l1 + unemp_l1 + realgdpgr_l2 + unemp_l2, data = macro_se)

summary(m_ardl)

coef(m_ardl)

theta1 <- m_ardl$coefficients[3]
phi0 <- m_ardl$coefficients[2]
phi1 <- m_ardl$coefficients[4]

est1 <- phi0
est2 <- phi0 + theta1*phi0 + phi1
est3 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1)
est4 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2)
est5 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2) + (theta1^4*phi0 + phi1*theta1^3)
est6 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2) + (theta1^4*phi0 + phi1*theta1^3) + (theta1^5*phi0 + phi1*theta1^4)
est7 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2) + (theta1^4*phi0 + phi1*theta1^3) + (theta1^5*phi0 + phi1*theta1^4) + (theta1^6*phi0 + phi1*theta1^5)
est8 <- phi0 + (theta1*phi0 + phi1) + (theta1^2*phi0 + phi1*theta1) + (theta1^3*phi0 + phi1*theta1^2) + (theta1^4*phi0 + phi1*theta1^3) + (theta1^5*phi0 + phi1*theta1^4) + (theta1^6*phi0 + phi1*theta1^5) + (theta1^7*phi0 + phi1*theta1^6)
estlr <- (phi0+phi1)/(1-m_ardl$coefficients[3])

se1 <- msm::deltamethod(~x2, mean = coef(m_ardl), cov = vcovHC(m_ardl))
se2 <- msm::deltamethod(~x2 + (x3*x2+x4), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se3 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se4 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3) + (x3^3*x2+x4*x3^2), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se5 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3) + (x3^3*x2+x4*x3^2) + (x3^4*x2+x4*x3^3), mean = coef(m_ardl), cov = vcovHC(m_ardl))
se6 <- msm::deltamethod(~x2 + (x3*x2+x4) + (x3^2*x2+x4*x3) + (x3^3*x2+x4*x3^2) + (x3^4*x2+x4*x3^3) + (x3^5*x2+x4*x3^4), mean = coef(m_ardl), cov = vcovHC(m_ardl))

Horizon <- c(1,2,3,4,5,6)
Estimate <- c(est1,est2,est3,est4,est5,est6)
se <- c(se1,se2,se3,se4,se5,se6)

impdat1 <- cbind(Horizon,Estimate,se)
impdat1 <- as.data.frame(impdat1)

impdat1$Estimate1 <- impdat1$Estimate
impdat1$Estimate2 <- impdat1$Estimate

impdat1$Estimate1[impdat1$Horizon>2] <- NA 
impdat1$Estimate2[impdat1$Horizon<=2] <- NA 

p3 <- ggplot(data=impdat1, aes(x=Horizon, y=Estimate1)) + geom_point(color="steelblue3") + 
  geom_line(aes(y=Estimate), color="steelblue3") +
  geom_ribbon(aes(ymin=Estimate-se*1.645, ymax=Estimate+se*1.645), 
              fill="steelblue3", alpha=0.2) +
  geom_hline(yintercept=0, linetype=2, linewidth=0.8, alpha=0.4) +
  ylab("Estimate") + scale_x_continuous(breaks = c(1,2,3,4,5,6),
                                        limits = c(1,6)) +
  scale_y_continuous(limits = c(-2,1)) + theme_pubclean() +
  theme(text = element_text(family="Times New Roman"),
        axis.line = element_line(color="grey", linewidth = 0.6),
        axis.ticks = element_line(color="grey"), legend.position = "bottom")


#### LP ####


m0 <- lm(ins_d0 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_se)
m1 <- lm(ins_d1 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_se)
m2 <- lm(ins_d2 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_se)
m3 <- lm(ins_d3 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_se)
m4 <- lm(ins_d4 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_se)
m5 <- lm(ins_d5 ~ SPENDING + insurance_spending_l1 + realgdpgr + unemp + realgdpgr_l1 + realgdpgr_l2 + unemp_l1 + unemp_l2, data=macro_se)

stargazer(m0,m1,m2,m3,m4,m5,type="text")

vcov0 <- coeftest(m0, vcov. = vcovHAC(m0, type="HC1"))
vcov1 <- coeftest(m1, vcov. = vcovHAC(m1, type="HC1"))
vcov2 <- coeftest(m2, vcov. = vcovHAC(m2, type="HC1"))
vcov3 <- coeftest(m3, vcov. = vcovHAC(m3, type="HC1"))
vcov4 <- coeftest(m4, vcov. = vcovHAC(m4, type="HC1"))
vcov5 <- coeftest(m5, vcov. = vcovHAC(m5, type="HC1"))

vcov0 <- vcov0[2,1:2]
vcov1 <- vcov1[2,1:2]
vcov2 <- vcov2[2,1:2]
vcov3 <- vcov3[2,1:2]
vcov4 <- vcov4[2,1:2]
vcov5 <- vcov5[2,1:2]

selpdf <- as.data.frame(rbind(vcov0,vcov1,vcov2,vcov3,vcov4,vcov5))

selpdf <- selpdf %>% dplyr::rename(SE = `Std. Error`)

selpdf$Horizon <- c(1:6)

p4 <- ggplot(data=selpdf, aes(x=Horizon, y=Estimate)) + geom_point(color="red3") +
  geom_line(aes(y=Estimate), color="red3") + geom_ribbon(aes(ymin=Estimate-SE*1.645,
                                                             ymax=Estimate+SE*1.645),
                                                         fill="red3", alpha=0.2) +
  geom_hline(yintercept=0, linetype=2, linewidth=0.8, alpha=0.4) +
  scale_y_continuous(limits = c(-2,1)) + scale_x_continuous(breaks = c(1,2,3,4,5,6)) +
  ylab("LP Cumulative Multiplier") + theme_pubclean() +
  theme(text = element_text(family="Times New Roman"),
        axis.line = element_line(color="grey", linewidth = 0.6),
        axis.ticks = element_line(color="grey"), legend.position = "bottom")

ggsave(plot = p1, filename = "IRF_INS_ARUK.JPEG", dpi = 800, height = 3, width = 3.5)
ggsave(plot = p2, filename = "IRF_INS_LPUK.JPEG", dpi = 800, height = 3, width = 3.5)
ggsave(plot = p3, filename = "IRF_INS_ARSE.JPEG", dpi = 800, height = 3, width = 3.5)
ggsave(plot = p4, filename = "IRF_INS_LPSE.JPEG", dpi = 800, height = 3, width = 3.5)


################################################################
